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Quantum states of a discrete breather are studied in two ways. One method in- 
volves numerical diagonalization of the Hamiltonian, the other uses the path integral 
to examine correlations in the eigenstates. In both cases only the central nonlinear- 
ity is retained. To reduce truncation effects in the numerical diagonalization, a basis 
is used that involves a quadratic local mode. A similar device is used in the path 
integral method for deducing localization. Both approaches lead to the conclusion 
that aside from quantum tunneling the quantized discrete breather is stable. 

PACS numbers: 31.15.Gy, 31.70.Hq, 05.45.-a, 63.20.Pw. 

I. INTRODUCTION 

Nonlinearity creates localized structures through a variety of mechanisms. In a discrete 
lattice, extreme displacement of one or a few atoms pushes them to high oscillatory frequen- 
cies outside the range of ordinary (linear) phonons, so that the energy sequestered in these 
large excitations does not spread throughout the lattice. Such an excitation, known as a dis- 
crete breather (DB) or intrinsic localized mode (ILM), is a kind of soliton, although different 
in many ways from the moving, continuum soliton of Russell ^ . The discrete breathers of 
primary interest to us, because of their relation to decay anomalies in doped alkali halides 
[21 , are not "intrinsic" in the sense that they are induced by excitation of the impurity, and 
translational invariance has been lost. For convenience, and to minimize acronyms, we refer 
to the excitations we study as "breathers." 

There is an extensive literature on these breathers and an introduction can be found in 
[^. Because of the nonlinearity, a good deal of the theoretical work in this field is numerical. 
For systems well-approximated classically this is not a problem, but treating the system 
quantum mechanically demands potentially serious approximations. 

An issue of experimental relevance is the lifetime of the quantum breather against decay. 
At the classical level, both analytical and numerical experience support an infinite lifetime, 
at least for the one-frequency breather j^. In a diatomic lattice there are two breather 
frequencies (in the phonon gap and above the optical phonons) and although there are 
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indications that combinations of the modes can allow energy to leave via the phonons, in 
our simulations this does not take place in any perceptible way. 

But classical stability says little about its quantum counterpart. Quantum tunneling 
bypasses classical constraints, and for a translationally invariant system, guarantees that 
bound states that are not of finite support become bands. Indeed such bands — as well as 
tunneling — have been studied for breathers in a homogeneous lattice |^,0,0|- Other breather 
quantization approaches have also been taken, including semiclassical and field-theoretical 0, 
|9[ . In the semiclassical approach one uses Einstein-Brillouin-Keller quantization on the KAM 
torus- confined classical trajectory. Naturally this yields a stable state, since, as for a 1-D 
state hidden behind a barrier, the classical trajectory does not sense the possibility of escape. 
Nevertheless, in principle it is entirely plausible that corrections to semiclassical behavior, 
which generally begin at order /i^, would include the possibility of decay (although in one- 
dimension the tunneling decay correction is the much smaller exp(— [positive constant]/^)). 

The focus of the present article is the stability of the quantized breather. At times it 
seems that the existence of even the classical breather is magical. Were this not a known 
phenomenon, one's first reaction on finding this numerically should be to debug the program! 
Why then should one expect this "magic" to operate in the quantum case as well? Perhaps 
to first order in h, where semiclassical considerations are a guide, but a priori there is no 
reason that quantum corrections should not lead to decay channels. 

Indeed, in a series of articles [13, U, [13, l3| it has been claimed that such decay takes 



place, and estimates for its rate in alkali halides give lifetimes of about 10 ns. Decay is found 
to occur irrespective of dimension and, as far as we can tell, the decay rate is of order h, 
surprisingly large. We recall the framework used in these articles, since to some extent we 
will take a similar approach. It is typical of the classical breather that a small number of 
atoms (often just one) vibrate strongly while the others hardly move at all. This strong 
vibration should not be treated by perturbation theory, and one replaces this dynamic, 
nonlinear object by something more tractable. Refs. jiol . 11, 0, replace it by a fixed. 



rapidly oscillating forcing term which couples — including nonlinear coupling — to the phonon 
field of the crystal. Two approximations are inherent in this: nonlinear interactions other 
than those with the central atom are neglected and the motion of the central atom is treated 
classically, subject neither to quantum laws nor to the back reaction of the phonon field. 

Recent work 3| takes a similar perspective but comes to different conclusions regarding 
decay rates. For the infinite lattice the authors find no decay at all, with departures from this 
situation for finite (lattice size) going like 1/A^. In any case, this is in strong disagreement 
with nanosecond-scale lifetimes in physical crystals. 

In the present article we treat the entire system quantum mechanically. We use two 
distinct approaches, having at their core fundamentally different approximations. In both 
cases however the unavoidably large impact of the breather is treated by introducing a local 
mode, either for comparison or for use as a zeroth-order system around which to perturb. 
As in other studies we drop nonlinear contributions except those that involve the particles 
most active in the breather. 

Our first method employs the Feynman path integral. An early triumph of that technique 
was the calculation of polaron properties |il5». ilQI]. and the principle lying behind its success 
was that almost all forces involved were linear. Using the path integral most degrees of 
freedom could be integrated, leaving the polaron with a self-coupling that was non-local in 
time and reflected the back-reaction of the lattice to its motion. All this was at the quantum 
level. By dropping all nonlinear terms except that of the most active breather atom we arrive 
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at a situation where the same method can be used. This also works if the active atom has 
only linear interactions, and allows us to compare properties of the breather to properties 
of a local mode, one whose motion is known to be fully localized. This is important because 
it is not enough to show the existence of eigenstates that correspond to the breather — 
one must show these eigenstates to be localized. What emerges from our calculations is 
that correlations of the breather-atom's motion with other atoms in the lattice drop in the 
same way as those of a local mode, as the strength of the appropriate coupling constant 
is increased. The path integral method has the additional advantage that its extension to 
higher dimension should be possible. (Note that our use of the path integral is unrelated 
to certain other uses of functional integration for solitons, for example as described in the 
work of Faddeev and Korepin 

The second method is direct diagonalization of the Hamiltonian. Essentially the same 
calculation was performed in j^, although they maintained translational invariance — they 
did not drop any nonlinear terms — and as a result obtained band structure. The difficulty 
of numerical diagonalization is that a cutoff on the size of the phonon Hilbert space is 
required. In a phonon basis, a quartic interaction couples to ever-higher phonon excitations, 
and perforce some of these must be neglected. Moreover, the cutoff must be fairly low. The 
elegance of the number-operator representation allows the analytic calculator to forget that 
the phonon spaces form a tensor product. Nevertheless, if n phonons are considered, with 
cutoff Hilbert spaces of dimension Nc, the Hamiltonian will live in a space of Nc" dimensions, 
potentially a computational disaster. Our innovation is to introduce a local mode into the 
zeroth order Hamiltonian, one whose frequency is roughly that of the breather in which we 
are interested. In this way, when going over to nonlinear coupling very few phonons are 
excited other than those of the local mode, and a low cutoff does little harm. 

Our conclusion from these calculations is that the quantum breather, at least in one 
dimension, is stable. (One-dimension, because of the intense distortion of the symmetry- 
breaking Jahn- Teller effect, is the case of physical interest to us 0.) From the path integral, 
the matching of correlations shows that the quantum breather has the same localization 
properties as a local mode, which is known to have an excitation profile (for atoms of the 
lattice) that drops off exponentially away from its source (as will be explicitly shown). For 
the direct diagonalization, the contribution of phonons other than those of the local mode 
to the true eigenstate is extremely small, on the order of 10"'^ for the parameters we use. 
We also address the issue of whether even very small contributions can lead to decay, which 
might seem to be implied by Fermi's Golden Rule. It turns out that coupling to a continuum 
does not always mean decay to the continuum, as was found in ^]. In the present article 
we give explicit examples and in particular examine the time evolution of a state initially in 
an eigenstate of the local mode — which is not an eigenstate of the full Hamiltonian. It does 
not decay. 

In Sec. mil we show that dropping nonlinear terms, except for those affecting the central 
atom, leaves the classical mechanics qualitatively unchanged. Following that, in Sec. II V[ we 
present the path integral approach. Direct diagonalization of the Hamiltonian appears in 
Sec. |V]and the associated time evolution in Sec. IVII A final section is devoted to discussion 
and conclusions. 
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II. DEFINITIONS AND NOTATION 

The system is a ring of + 1 unit-mass atoms with Hamihonian 

N 

H = Y.{^Pk+'2^s4 + i^lixk-Xk+i)^ + \Xxl} , xeR, (1) 

A;=0 

Periodicity is expressed through Xq = xtv+i and mod-(A^+l) addition for atom labels. This 
is the nonlinearity studied by jsf. Another system, closer to our own models, is 

N 

H = ^{yi + luj^xl + luj^ixk - Xk+if + jX{xk - Xfc+i)^} , (2) 

k=0 

also periodic. To avoid being awash in subscripts and superscripts we distinguish these by 
context, rather than by labels attached to H. As remarked, when a breather is present, one 
typically has a single atom where most of the energy is concentrated, with the other atoms 
relatively still. We take the dynamic atom to be #0 and consider an alternative Hamilto- 
nian in which all nonlinear terms that do not involve Xq are dropped. The corresponding 
Hamiltonians are 

N 

H = Y.{H + i^'^' + l^oi^k - Xk+i f] + \\xt , (3) 

fc=0 
N 

^ = I] (bfc + l^s4 + Woi^k - Xk+iY] + IX [{xo - x,y + {xo - xn)'] . (4) 

k=0 

For each of these we introduce a local mode Hamiltonian, namely a linear system which, at 
a level to be explored below, behaves similarly to those in Eqs. Q and They are 

N 

^0 = ^{IpI + + l^lixk - Xk+iY} + \ujIxI , (5) 

k=0 

N 

Ho = {|Pfe + l^s4 + l^^U^k - Xk+if} + ^ujI [{xo - xif + {xo - xnY] ■ (6) 

k=0 

As the zero subscript suggests, these will be the zero-order perturbation Hamiltonians in 
our numerical diagonalization. The perturbations will then have the form 

Vj = \\xl - \u:\xl , (7) 
Vi = iA [(a;o - Xi)^ + {xq - xn)'^] - ^uf [{xq - xif + {xq - x^f] , (8) 

for the Hq's of Eqs. ^ and (jHI), respectively. 

All these systems have reflection symmetry about atom-0. Define P by 

P : A; ^ (iV - + 1) , A; = 0, . . . , iV_ with iV_ = [^] , (9) 

and the square brackets represent "integer part of." Classically an initial condition having 
this symmetry will retain it. Quantum mechanically, since Vj shares this symmetry, only 
phonons even under P need be considered. 



5 




500 



^tom 



30 



FIG. 1: Kinetic energy as a function of time. Only one figure is shown, as the others are indistin- 
guishable at this level of display. 

III. EFFECT OF DROPPING NONLINEARITY FOR "NON-BREATHER" 

ATOMS 

Since we drop nonlinear terms that do not include Xq, we first study the impact of that 
change on the classical mechanics. The dynamical systems to be compared are 

• of Eq. (P), "All nonlinear." 

• H oi Eq. (j2)), "One nonlinear" {xq in the equation). 

• Hq of Eq. (jni), "Local mode" (around xq). 

We also make similar comparisons for nonlinear coupling: 

• of Eq. (0), "All nonlinear." 

• H oi Eq. (0}, "One nonlinear" {xq couples nonlinearly to its neighbors). 

For the self-coupling cases, we first exhibit the kinetic energy as a function of time for 
a ring of 30 atoms. As for our subsequent graphs we have shifted the central atom to the 
middle of the ring to avoid cutting the breather. The kinetic energy is a clear way to see 
the breather and is shown in Fig. ^ In Fig. El we present position as a function of time 
for all three cases. Finally to illustrate the fact that differences do remain, we show the 
Fourier transform of the position of the central breather particle. Fig. El shows intensity, 
and it is clear that when all atoms see the nonlinear force there are more pronounced higher 
harmonics in the motion of the central atom as well. All show activity in the phonon band, 
notwithstanding the fact that all are stable. 
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FIG. 2: Position as a function of time. The motion of each atom is displaced by one, with the 
active breather atom in position 16. The figures ihustrate symmetry and rapid dropoff as well 
similarity among the various dynamical systems, a) Full nonlinear system, b) Only nonlinearity 
in breather atom, c) Local mode at breather atom. 
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FIG. 3: Fourier transform of the position of first atom. The "One nonlinear" plot is displaced 
downward by 3 and the "Local mode" plot by 6. Dotted lines indicate the span of the phonon 
band [which does not begin at because of the ujs of Eqs. through ©]. 



The reason for the relative indifference of breather behavior to nonlinear couplings not 
involving #0 is that in any case atoms other than 7^0 have small amplitude and hardly 
feel the nonlinearity. It is certainly true that for smaller A, when the breather spreads over 
a large number of atoms (as occurs in several cases in |l9|), this approximation would be 
inadequate. However, in this article we consider only the highly localized case so as to 
focus on the issue of in-principle quantum stability. Fig. EJ based on the full Hamiltonian 
Eq. (P), shows how small the nonlinear contributions to the energy are for any but the 
central breather atom. What are displayed are the time averages of lolx^ and . It is clear 
that even for the nearest neighbors of the highly oscillatory atom the nonlinear contribution 
to the energy is negligible. 

For the case of the truncation of nonlinear coupling, the parallel is not so close. The 
Hamiltonian, Eq. (j21), gives rise to interesting patterns of behavior not seen in the self- 



coupling case. There is first the phenomenon of moving breathers |20l . |21L l22l |23L |2 
2^, which obviously will not occur when only one site and its neighbors are subject to 
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FIG. 4: Comparison of linear and nonlinear contributions to the energy. Nonlinearity is retained 
for all atoms. The upper figure shows the time averages of iOgx"^ and Xx^ for all atoms. The lower 
figure shows a subset of the same information, but on a scale large enough to see that even for the 
nearest neighbors of the most active atom the nonlinear contribution is extremely small. 
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FIG. 5: As in Fig. [21 The figures illustrate symmetry and rapid dropoff as well as the differences 
between the full and truncated Hamiltonian. a) Full nonlinear system, b) Only nonlinearity in 
breather atom coupling. 



nonlinearity. Movement can be suppressed by taking appropriate initial conditions (equal 
and opposite position displacements) so that the stationary breather consists of a pair of 
atoms oscillating with respect to each other. To reproduce this in a truncated Hamiltonian 
we would have to keep the nonlinearity for at least two atoms, which, for our later numerical 
diagonalization, puts a strain on computer memory. We therefore stayed with the form 
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Time -> 



FIG. 6: Positions as a function of time for a "wanderer." The motion is governed by the same 
equations that yield Fig. |SJ but the initial conditions are different. Here one atom is displaced by 
a relatively large amount. 



which provides a local, nonlinear structure. The logic of [10|, |11|, ll3[ predicts decay of 
this system when quantized, so that finding the eigenvectors of this system to be localized 
provides a further counterexample to that logic. In Fig. El we show the contrasting behavior 
of position as a function of time for the true breather and the truncated one. Apropos moving 
breathers, in Fig. IHl we present an illustration of a typical example, wandering around the 
ring in a seemingly random way. 



A. Quantum perspective 

Quantum mechanically, dropping nonlinearity away from the center of the breather has 
a profound effect. It represents a loss of translational invariance, and, because of quantum 
tunneling, a qualitative change in the spectrum of the system. It follows that, in principle 
at least, for a translationally invariant system all breathers are nonlocal. The associated 
band structure for breather levels was found in ^ and Ref. [i^ has an extensive discussion 
of the significance of translational invariance. 

It is therefore important to state what we mean in asserting that the breather is localized, 
and why our assertion is physically significant. If one considers a deeply bound atomic level, 
say the innermost shell of Cu when the atom is part of a crystal, that level — in principle — is 
part of a band. Practically, however, the band is irrelevant and the measurable properties of 
this level are independent of the near-zero tunneling probabilities that induce the theoretical 
band structure. In our work, we will consider breathers under circumstances where they are 
far from classical instability, so that quantum processes that connect to distinct classical 
breather states (e.g., the breather translated by one atom) would be of extremely small 
amplitude So we are saying that except for this small tunneling amplitude the breather 
is stable. Physically this means that the breather can have an extremely long lifetime, 
contrary to the claims made in the literature that we have cited U, [13, We 
emphasize that in that literature translational invariance is also dropped. 



A second consideration is that the physical systems of particular interest to us [2] do not 
have translational invariance. The breather formed by Jahn- Teller distortion at an impurity 
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FIG. 7: Selected atoms on the ring. For convenience the notation q is used for xq. 



certainly has a preferred origin, and it is stability against low order phonon decay that 
concerns us. 



IV. THE PATH INTEGRAL APPROACH 

Feynman's use of the path integral for the polaron ^ exploited what has come to be 
one of the most useful features of the path integral: the propagator for a quadratic degree 
of freedom coupled linearly to something else can be evaluated explicitly, leaving only a 
self-coupling of that "something else." The penalty is the nonlocality in time of the self 
coupling, the reward the reduction of the problem to a single degree of freedom. 



A. Setting up the path integral 

As discussed in Sec. IllH for appropriate parameters, all atoms but one have negligible 
nonlinear energy contributions, and we drop the "Ax^" for all but one of them. This allows 
us to integrate all but the single atom, around which the breather is centered. As in Eq. Q, 
we take that atom to be #0. As a result we consider the Lagrangian 

^ N+l ^ N+l ^ N+l ^ 

C= - ^Xl- -Ujl^iXn- Xn+lf - -U^ XI ^« ~ 4^0 + f^M^m + ^N+l-m) ■ (10) 
n=l n=l n=l 

We have introduced an additional term, fiXo{xm + XN+i-m)- It is easy to carry along and 
will later allow us to study correlations. The derivative (at yU = 0) of the (appropriate form 
of the) propagator provides a (O-m)-correlation. We will take particle-m to be distant from 
0, typically about 1/3 of the way around the ring. The use of a pair [m and (A^ + 1 — 
m)] maintains mirror symmetry ("P") with respect to particle-0. We will also distinguish 
between the ring, all + 1 atoms, and the chain, atoms 1 through A^. See Fig. [7| 

The Lagrangian C can be written in the form C = Cq+Cc+Ci. Here Cq is the Lagrangian 
of particle-0 

^0 = \e - U + \A - \q\ (11) 
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where we have adopted the notation q for xq. Cc is the Lagrangian of the chain 

1^1^ 1 ^ 

n=l m,n=l n=l 

with Jmn = '^Smn — Sm+i,n — ^m.n+i ("^, u = 1 . . . N) , the Jacoby matrix. The last term £/ 
is the interaction between the particle-0 and the others 

Cl = UqQ (Xi + Xn) + M{Xm + XN+l-m) ■ (13) 

We next perform the path integral for the chain and arrive at an effective action and 
path integral for q alone |l5l Eg^ .29. . .3Q] . Normal coordinates for the chain are defined by 
Q = Sx, where Smn = \/2/ {N + 1) sin [vrmri/ {N + 1)] is a unitary matrix, and {xi, . . . , xn) 
is considered a column vector. This leads to 



1 ^ 

+ = 2 E - ^IQI + fn{t, q)Qn} • (M) 



n=l 



Here {f2n} are the chain normal-mode frequencies 



nl = + Aul sin^ j^P^^ , (15) 



and /„(t, q) = qr^ [1 - (-1)"], where 



T„ = \ ciJn sm h u sm . (lb) 

The path integrals over the normal chain coordinates Qn can easily be evaluated, so that 
one has a propagator that is a function of the endpoints of all degrees of freedom. It is 
convenient to take the matrix element of this object in the chain ground state and divide 
by corresponding matrix elements of the unforced chain oscillators. This yields an effective 
action for the chain 



/N 1 „ Af 

L n=l J n=l 

where = J C^, 'j = C,I, and (e.g.) ^pQ is shorthand for Hn ^o^"'*(Qn)' ^^"^ '^o^\Q) = 
(to /nh)^^* exp(—ujQ'^ /2h) is the ground state of a harmonic oscillator. The chain oscillator 
modes are subject to a forcing term due to the so-far-unintegrated q(t). 

The result of these standard calculations is that the propagator for particle-0 takes the 
form 

g {q^, T- q\ 0) = y Pgei(^o+^^«) , (18) 
where the action 5*0 arises from the original L sans chain terms, specifically 

^0= /^4^^'-(^o + ^-.')^'-I^4- (19) 
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The effective action, S^g, is the result of the integration just described over chain degrees of 
freedom. It is given by 



S^^= [ I dtdsK{\t-s\)q{t)q{s), (20) 
Jo Jo 

with 

^ ^ " ^^sin(^) ^ ^ 



n=l,3,. 



We shall study the (matrix element of the) propagator in the stationary phase approxi- 
mation; that is, we focus on the action (including S^n) along the extremal "classical paths." 
The spectral expansion of the propagator [as a function of all variables, i.e., before the 
chain-variable integrations and the division of Eq. (fT7j) ] is 

G{xl x'l, ...,t;x'„ x[, ...) = J2 ^1' • • •) eM-^Eat)Ki4, x[, . . .) , (22) 

a 

where a is an eigenstate label. This implies that the operations of Eq. ()17|) (the integrations, 
and the division, which leads to energy-phase corrections) lead to 

g {q", t- q\ 0) = exp {+i ^ fi„t/2) ^ 0„(g") exp(-zE,t)C(g') , (23) 

a 

where 

0„(g) = j ^o{xi, . . .y^aiq, xi,...)dxi... , (24) 



and \l/o is the ground state of the chain. Note that Ea in Eq. ()23|) is the total energy of the 
ring, while Eq = X] is the total ground state energy of the chain. Going to imaginary 
time, t = —iT, Eq. becomes 

g {q", T- q\ 0) = J] 0„(g") exp [-{E^ - Ec)T] MqT (25) 

a 

This is the starting point for many calculations. The simplest would be to get the ground 
state energy of the chain by letting T —>■ oo. For numerical work this turns out to be delicate. 
What one would do is find a semiclassical expression for the total action in Eq. (jl8|) and 
consider its large T limit. Unfortunately S tends to a constant, and the energy emerges 
from its derivatives. The same occurs for the simple harmonic oscillator. Thus for Csho = 
mx^/2 — mi/^x^/2, the propagator from a to 6 in imaginary time T is 



1 d'^S 

GsUb,T;a) = \/^^exp(-^. 



2ti sinh uT 



For large T the dominant term in the exponential approaches a constant, so that the energy 
comes from the prefactor, a/ d'^S/dadb, which as the explicit form above shows is essentially 
1/Vsinh z/T, yielding the correct ground state energy, u/2. This means that one needs 
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FIG. 8: The classical local modes are U£{n), with i the mode label and n the atom. Shown in the 
figure is the base-10 logarithm of the absolute value of mode-0, the local mode. Parameters are 
ioi = 2 and ujq = 1 [cf. Eq. ©]. Clearly there is an exponential dropoff of excitation for the atoms. 

extremely accurate values in a second derivative of the action, exponentially small in T. For 
the action itself this precision was attainable in all cases, and for a single harmonic oscillator 
the second derivative was also reliable, but for the full ring we did not find that the second 
derivative had the needed precision. 

In any case, evaluating the energy is not our primary goal. We know that the Hamiltonian 
has eigenstates, but we do not know whether they resemble breathers. We also expect that 
some states, not necessarily eigenstates, but at the least metastable states, should resemble 
the classical breather state. This much is guaranteed by semiclassical considerations. Our 
goal then is to show that the actual eigenstates resemble the breathers. 

To see how to do this, we will consider first the local mode, a quantum state that we 
know to be localized. Specifically, the classical mode function has a sharp spatial dropoff 
and the quantum state represents atomic oscillations with that pattern. For reference we 
show the classical local mode for ui = 2 and cuq = 1 in Fig. |H1 [a solution of Eq. ©]. (The 
value of Us does not affect the function.) 

Our goal is to ascertain properties of (pa and Ea from the semiclassical approximation to 
G. This goal will be attained in Sec. IIVCI However, before actually doing this calculation we 
discuss methodology. Both for the linear and nonlinear cases, implementing the semiclassical 
approximation presents technical challenges, which we now discuss. We will also indicate the 
extent to which independent confirmations of our technique are available, taking advantage 
of known asymptotic properties. 

B. Methodology 

The first observation to be made is that the development from Eq. ()10|) to Eq. ()25|) goes 
through in exactly the same manner, with the same K, etc., for the quadratic local mode 
as for the nonlinear excitation, with the sole exception of the final form of 5*0, which, for 
the local mode case, has a term u)fq'^/2 instead of Going over to imaginary time, our 

problem (both quantum and classical) revolves about the following quantities: 

^to.a, = ^ dt!^^e + ^{u;'^+2u;l)q' + U{q)\- ^ dt ds K{\t - s\)q{t)q{s) , (27) 



13 



where U{q) = a;^g^/2 for the local mode problem. The same form holds for the nonlinear 
problem as well, but with U{q) = S'^otai is the total classical action, and in the path 

integral appears as exp(— S'totai) (note the minus sign, a result of the t — > —it transformation). 
The self-interaction kernel becomes 

The classical equation of motion follows from the usual variational methods and is 

q - {^l + 2ul)q - — + 2 K(\t - s\)q{s) = . (29) 



For the semiclassical approximation one must solve Eq. ()29p . We used different methods 
for the linear and nonlinear cases. We discuss them separately. 



1. Linear, non-local in time propagators 

For U{q) = ujfq'^/2, Eq. can be discretized, say with = ke, k = 0,1, . . . , M + 1, 
e = T/M. q becomes a vector, q^ = q(tk) {k = 1, . . . , M) and the entire equation has the 
form Bq = qo, with B a matrix consisting of three parts: the second derivative operator, 
2/e^ on the diagonal, — above and below; a diagonal matrix proportional to {Ug + 
2u!q) and a matrix (which by a slight abuse of notation we call) K whose substantial non- 
diagonal components reflect the nonlocality in time. From Eq. (j29p it may not be evident 
how a nonzero right- hand- side (go) comes into the picture. It arises from the boundary 
conditions. When discretizing, the first and last rows of the second derivative operator call 
on g-components outside the range 1, . . . , M. These other components are the boundary 
values (call them a = ginitiai and b = qan^i), so that go(l) = — a/e^ and qo{M) = —b/e^ 
(other components of go are zero). This method of dealing with the two-time boundary 
value problem is a simpler version of what is done in 3l|. It follows that the solution of the 



nonlocal-in-time equations of motion is g = B~^qo. It is also immediate that for the linear 
equation the action along this "classical path" is given by S* = [q(t)q{t)/2]\^ , although one 
can also evaluate S by exphcit integration. 

The substantial nonlocality in time means that the matrix K is not sparse, putting 
a limit on the smallness of e, hence on the accuracy of the action. This difficulty was 
overcome by the following device. For each such e, call S{e) the action that results from 
the discretization and matrix inversion just described. Considered as a function of e, one 
can write S{e) = S{0) + e5"(0) + e^5"'(0)/2 + . . ., where the prime indicates derivative with 
respect to e. By evaluating ^(e) for several values of e one can then extrapolate to zero. For 
example, the simplest such extrapolation, for two values, ei and €2 (and associated 5*1 and 
5*2), gives 5(0) = (€25*1 — €182)/ {€2 — ei). We typically used three values, chosen as small as 
possible, but far enough apart to keep denominators from being too small. As we indicate 
below, this gave excellent results when tested against quantities that could be independently 
calculated. The following tests were performed: 

• Simple harmonic oscillator. For the energy this requires taking the second derivative 
with respect to position. This was done by evaluating G, hence S, at several points 
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and using 



dadb 



S{a + 6,a + 6) + S{a — 6,a — 6) — S{a + 6,a — 6) — S{a — 6,a + 6) 



(30) 

For this formula the extrapolation to 5 = was also used (as for e above), which is 
to say, several values of 6 were used. Furthermore, there was also an extrapolation 
to T = oo (by fitting the logarithm of d'^S/dadb as a function of 1/T). Results were 
accurate to better than one part in 10^. 

Asymptotics of 4>o{<l) [defined in Eq. pij) ]. This can be calculated by non-path integral 
methods as follows. All eigenstates — both ring and chain — are Gaussians. For the 
chain the frequencies are nt^'''^ = + Auism\n7r/2{N + 1)), n = 1, . . . , A^. For 
the ring, for the case ui = 0, they are Qn""^^ = ^y + Au^ sin^(n7r/(A^ + 1)), n = 
0, . . . ,N. If a;i 7^ 0, they are easily evaluated numerically. Aside from normalization, 
00 is therefore given by [cf. Eq. (j^ ] 

0o(g)= / n ^Qi"""^ (-Efc^l"""^Qi"^'"^V2) exp (-E,^^^^^gJ^^^^^^ . (31) 

^ n=l 

Given the transformations from the original variables, {q,Xi, . . . ,Xn) to the chain 
and ring variables (known analytically for the chain, numerically for the ring), it is 
straightforward to get the coefficient of in 0o- Let the matrices ?7*-""*^^ and f/(<^'^=''") be 
defined by Q^""'^'' = U^^^^'^Xi and Q^'"'"''^ = Yli ^u'^'^'^^t (with {Q} the appropriate 
normal coordinates). The f/'s can be taken to be real. Let 



N 



Uk 



uj^^, n^^J2^tM, (32) 



N N 

9n ^ (33) 

A;=0 e=l 

N N N 

W^„. ^ S^^^Q^^-^ + E E E ut^^'^^ut^^Qt^uj^'^ujC'^ ■ (34) 

A;=o i=i e'=i 

(Note that t/C'^hain)! ^ f/(chain) rj.^^^^ coefficicut of in the exponent of (pl (call 
it C) is 

( = Qu-9^W-'9. (35) 

To evaluate ( by the path integral method, the action was evaluated with various 
endpoints, a. Since —S{a,T; a) (the action from a to a in time T) appears in the 
exponent of the propagator, if the semiclassical method and our extrapolations are 
valid, S should be quadratic in a. This was tested by fitting \/S against a. We did not 
do this for a variety of T values since a single large T was sufficient. Large in this case 
means with respect to the smallest spectral gap, and can be discerned numerically by 
a stabilizing of the slope in the fit. The desired result is that first, the fit of against 
a should be a straight line, and second, that the square of the slope of the line should 
agree with (. We will not reproduce the graph of VS against a since there is little to 
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note beyond its being an excellent fit. The values of slope- squared and ( agreed to 
better than 1% for min(ti;i, 0;^) > 1 and by no more than 3% when these quantities 
became smaller (they impact the spectral gap, which can force T values too large for 
good accuracy). 



2. Non-linear, non-local in time propagators 



With a nonlinear term in the action, the classical nonlocal two-time boundary value 
problem cannot be solved by matrix inversion. The equation to be solved is 



q - {ui + 2uj',)q - Xq' + 2 K{\t - s\)q{s) = . 

Jo 



(36) 



Our method is an extension of Feynman's approach [1^ in which the introduction of an 
auxiliary variable eliminates the nonlocality. Then one can use standard numerical methods 
for solving the two-time (even-if-nonlinear) boundary value problem. Recalling the definition 
of K, Eq. (j^ . one can define auxiliary variables Zn by 



^ coshfi„ (f - |t - s\) 
as — _ rr.^ — Qys) ■ 







fi.sinh(^) 



Eq. (jHHjl becomes 



q - [uj'. 



+ 2ul)q- \q^ + 2Y,r^Zn = 0. 



Taking two derivatives of Zn leads to 



nizn - 2q{t) 



n = 1,3, ... ,2 



N-1 



+ 1. 



(37) 



(38) 



(39) 



One non-local equation has been replaced by a larger number of local ones. There is but a 
single pair of boundary conditions: g(0) = a, q{T) = b. The conditions on z are forced by 
self-consistency arising from Eq. (j?fj) . which incidentally also imply z„(0) = Zn{T), for all 
n. It's amusing that with the quartic interaction replaced by a quadratic one, this linear 
system is equivalent to the set of classical linear equations that apply to the ring [i^, with 
the boundary conditions inherited from the ground state averaging [which is of course where 
Eq. ()36|1 (made linear) comes from in the first place]. 

For numerical solution of Eq. (j36p we applied a variation of this method. As a function of 
its argument, K{u), although a sum of many hyperbolic cosines, actually (for large and 
T) bears a strong resemblance to a single such function. For the parameter range of interest, 
one can, by judicious choice of cUed, bring the difference between K and its approximation, 
Ki{u) = K{0) cosh{u^gu), to about 10~^ (integral of square of difference). With a sum of 
two hyperbolic cosines (and 3 adjustable parameters) one can do much better, but as we 
shall see below a single hyperbolic cosine was sufficiently accurate for excellent minimization 
of the action, i.e., solution of the classical motion. 

With a single hyperbolic cosine in the nonlocal portion of the equation, only a single z 
need be defined [cf. the derivation of Eq. (jSHI)] so that our nonlocal equation becomes a pair of 
second order ODE's, with boundary values for one of them and a self-consistency condition 
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for the boundaries of the other. Using the optimal uj^g defined in the last paragraph, the 
equations are 

q = {ul + 2ul)q + \q''-2z, (40) 



T 



coshcj^fj — \t — s|) 



z = uji^z - 2q{t) , with z{t) ^ / ds l-\^^^Lq(s) . (41) 



a;.ffSinh(5iif^) 

The boundary values of q{t) are given. For z one proceeds iteratively. For given z{0) (which 
automatically equals its value at T), the boundary value problem can solved using MATLAB's 
program "bvp4c" 3^1 . Using the associated solution, q(t), one can recompute z{0). In effect 



one has a function mapping z{0) into a new value, and again, standard numerical search 
techniques can be use to find a fixed point of this map. That fixed point then provides a 
solution of the original boundary value problem for q. 

In principle with a better approximation for K, using say two hyperbolic cosines, the 
parameter space of boundary value mappings becomes 2- dimensional. This, however, was 
not the way we proceeded. 

Since the true solution of the classical nonlinear, nonlocal problem minimizes the action, 
5*, it is possible to improve the solution by modifying the q derived from the process given 
above in such a way as to reduce S. The class of functions to add to q{t) for this purpose 
can be narrowed by the following consideration. The potential that we study is an inverted 
oscillator, linear or nonlinear. Therefore if the time interval for going between two not- 
too-small boundary values is large, for most of that time interval the particle will be near 
zero, with zero velocity: the path is thus an instanton and, except near the endpoints, will 
be exponentially small (in T). Profitable, i.e., S'-reducing, variations of q will thus have 
the same shape. Our basic variation consists of adding and subtracting hyperbolic cosine 
functions with varying amplitudes and angular frequencies. 

Thus, starting with the approximate- if solution to the nonlinear boundary value problem, 
we allowed modifications of the sort just described. The changes in S that resulted were 
extremely small. The first such correction was on the order of 10~^ of the action, and 
subsequent reductions were of order 10~^^. No other functional forms (e.g., multiplying 
instanton-shaped curves by oscillatory functions) gave any improvement at all. The results 
presented below all use this method of optimizing 5*. 

It is also possible to perform some of the checks made on the linear problem. In par- 
ticular the asymptotic form of the wave function for a quartic anharmonic oscillator is 
ip ~ exp (—[positive const ■ q^]). With various A we found the action for large boundary 
values of q. There was a very good fit to the cubic. In Fig. |Hlwe plot various powers of the 
rescaled action, showing that the dropoff for what we have called (poi^l) is close to cubic (the 
exact wave function is also not exactly cubic). 



C. Localization 

The local mode represents a localized excitation. Fig. |H1 shows the classical mode oscil- 
lation amplitude, hence the shape of the phonon; in particular it indicates an exponential 
dropoff with distance from the center of the mode. We will use our path integral formalism 
to show how it too reflects the localization of the "local mode." Then we will use the same 
technique to establish the localization of the nonlinear breather wave function. 
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2 4 6 8 10 

boundary value ("a") 



FIG. 9: Action as a function of boundary value. Since the wave function is essentially exp(— S"), 
the graph gives the its spatial dropoff. Like the wave function for a single quartic oscillator, the 
behavior of (poiq) is essentially exp(— g^). For this figure, A = 10 and T = 3. 



Recall the fictitious coupling, fixo^Xm + XN+i^rn), inserted in the Lagrangian of Eq. (fTUj). 
Site-m is far from site-0, where the large oscillations of the breather are taking place. Com- 
paring the action for small fi and the action for zero fi provides a correlation function in 
the following way. The imaginary-time version of Eq. (fH?)) is 



g{qf,T;q\0) = j Vqe 



(42) 



where there are slight differences from the definitions given in Eqs. (fT^ through (j?T|) . Specif- 
ically, the action (S'totai) is now given by Eq. (jTfj) . with K as defined there [note that Eq. (j2Z|) 
includes both the cases of quadratic and quartic potentials] . We emphasize that the fictitious 
coupling appears in the definition of r„ [Eq. flTH|)]. and only there. 

Thinking back to the propagator before the integration over chain mode ground states, 
we consider the derivative dG{a, —iT; a) /dfi\^=Q. The importance of this derivative arises 
from the relation 



J exp(— S* + fj,A)dx I J exp{—S)dx 



(43) 



where exp(— 5) is a weight for averaging. Since we are taking d/dfi after int egra ting, we are 
getting an average of the (ga:OT)-correlation in the chain mode ground state jSj]. 

Note too that because of our use of the semiclassical approximation (which happens to 
be exact in the linear case), study of dG{a, —iT; a)/dii\^=Q is essentially the same as study 
of a{a,T,g) = dS{a,T)/dn\^=o, where "5"' is the imaginary time action, a is the common 
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initial and final endpoint, and T the imaginary time. The quantity g, implicit in S{a,T), 
parameterizes the relevant particle-0 enhancement; for the linear local mode it is Ui and for 
the quartic case it is A. 

Consider first cr(a, T, ui) = dS{a, T)/dfi\^=o for moderate T (such that states other than 
the ground state survive in the spectral sum for Q), for small uJi and for large a. Because 
a is large, the important terms in the spectral sum will not be those of lowest energy, but 
those that permit large excursions of q [cf. Eqs. fE!?jl and (j^ ]. However, since ui is assumed 
small, there will be no local mode. So pulling q far from its equilibrium position, pulls all 
atoms far from their positions, and the correlation of atom-0 (g) and atom-m {xm) should 
be large. On the other hand, suppose uJi to be large (with a still large and T moderate). In 
that case there is a pronounced local mode and pulling q away from equilibrium has little 
impact on Xm- For the local mode, q can have large excursions while other atoms hardly 
move. Therefore one expects a to be small. In Fig. ITIIk . the lowest curve shows just this 
behavior. The boundary value of q, a, is 4 and for small uJi (no local mode) there is a large 
correlation with the motion of atom-m (0 is at the top of the figure). As uji increases, this 
correlation shrinks. 

By contrast, if a is small, the variation of ui has little effect. This can be understood 
as follows. The requirement on the endpoints of q imposes little demand on any other 
coordinates whether or not there is a local mode. 

We next turn to the nonlinear case. Here we do not have a priori knowledge of the wave 
function but can use the correlation function as a test of localization. The behavior of the 
correlation function, as a function of a and A exactly parallels that of the linear local mode. 
This is shown in Fig. ITUb . For small A and large a there is no breather and, as for the 
small Ui case, the demand for a large excursion of q forces a large excursion of Xm- And 
now the central observation: for large A, forcing q to be large has almost no impact on 
Xm, exactly as for the quadratic local mode, from which we deduce the localization of the 
breather excitation. We remark that corresponding values of ui and A are related by A ~ oof. 
Finally we mention that from Fig. ITUb alone it is difficult to tell whether dS/ d/j, is tending to 
a constant or to zero. In Fig. ITTlis a plot of {dS/d^Y versus 1/A, in which the extrapolated 
value is close to zero. 

Another check of correlation function behavior is its time dependence. For a localized 
stable state and fixed a (the boundary value for q) dS/dfi should decrease with time. This 
follows from consideration of the propagator given in either Eq. (j22|) or Eq. (j23j) . For large 
T and a there is a competition between the terms in the expansion, with large T favoring 
states of lower energy and large a favoring states that drop off most slowly in their spatial 
coordinate, which in general will be of higher energy. As T increases lower energy states 
are increasingly important in the mix and the correlation correspondingly reduced (note 
that for the correlation the overall magnitude of the propagator, and in particular the factor 
exp{—EQT), drops out). In Fig. [T21 this decline can be seen for both the quadratic local 
mode and the quartic breather. On the other hand, if the quantum breather were decaying 
in time the correlation should increase, since the true ground state of the system (if the 
breather were not stable) would resemble the A = case, for which the correlation is high. 

V. DIRECT DIAGONALIZATION 

The full Hamiltonian governing our system is given by Eq. Q [or Eq. (@))] and it is natural 
to seek the properties of the system by numerical diagonalization of H in an appropriate 
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FIG. 10: dS/d^ (essentially a correlation function) has markedly different behavior for large and 
small endpoint ("a") values for q. Both for a linear local mode (as a function of loi) and for the 
breather (function of A), a large value of a demands relatively large correlations when neither a 
local mode nor a breather is present (wi ~ A ~ 0), but that correlation is wiped out for large 
parameter (A or cji), for which the breather or local mode is effectively decoupled from the rest of 
the ring. 
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FIG. 11: Plot of {dS/dnf versus 1/A for a = 4. Same data as Fig.HHt). 



basis. Indeed this is the approach of [5|. The difficulty lies in the fact that a cutoff in 
phonon number is necessary in order to keep the problem finite. Since the overall Hilbert 
space is a product of individual phonon Hilbert spaces the total dimension is the product 
of that of the individual subspaces, so that the size of the matrix representing H is large. 
Moreover, even when one has the eigenstates of H (which by definition are stable), one still 
needs to show them to be localized in order to identify them with breathers. (Evidence for 
localization is also given in [^.) 
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FIG. 12: dS/dfi as a function of time. The figures show a dechne with fixed a and couphng 
parameter (ivi or A) as T increases. This is due to the increasing weight, with increasing T, of the 
ground state in the propagator. Cf. Eq. 



To reduce the burden of dimensional proliferation we used two devices. First, we introduce 
a second system for comparison, namely the Hamiltonian Eq. (jSJ [or Eq. ©], representing a 
local mode. The idea is that if uji (the fictitious self-coupling at site-0) is chosen appropriately 
the frequency and to some extent the shape of the breather (induced by a quartic with 
coefficient A) can be reasonably approximated. In that way, the perturbation, has minimal 
effect, and only a small number of non-localized phonon states enter the eigenstate associated 
with the breather. See Fig. El for an illustration of the efficacy of this method. The 
second device takes advantage of the reflection symmetry embodied in P [Eq. ©]. Since 
the perturbation {Xxq/A — ulxl/2) commutes with P we can focus on states in the same 
symmetry class as the local mode, namely those with P-eigenvalue 1. This reduces the 
number of phonons by almost a factor two, lowering the Hilbert space dimension to a bit 
more than the square root of what we would otherwise need to consider. Note that using 
eigenstates of P means that we are not using "traveling waves" for the phonons, and indeed 
one should no longer expect these to be the preferred basis once translational invariance has 
been dropped. 

Some details of the calculation play a role in interpreting the results and we present them 
here. The normal modes of our monatomic ring [Eq. (0)] satisfy the classical equation of 
motion, 

Xk + uj^Xk + ujI {2xk - Xk-i - Xk+i) + Sko^^lxk = , k = 0, . . . ,N + 1, (44) 

with mod-(A^ + 1) addition and 6jk the Kronecker delta. For lvi = 0, this is trivially solvable, 
while for non-zero Ui there is a local mode. Corresponding equations hold for the nearest- 
neighbor nonlinear-coupling model [Eq. ©], but since the principles are the same we do not 
present the equations in detail. With uq = 1 and Us of order unity, a typical spectrum is 
close to that of the ui = case, except for a single mode with frequency high above all the 
others. An example is shown in Fig. El 

At the quantum level the Hamiltonian is written using the classical modes. Call 
the frequency of the mode, u^^^ the coordinates of the mode [i.e., u^^^ exp{iil£t) solves 
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Local mode as a zero-order Hamiltonian 
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FIG. 13: Advantage of using the local mode for perturbation theory. Let Hq = + x'^/2, 
V4 = V2 = uJix'^/2, minus 2 = V4 — V2. For this graph we use A = 8 and uji = 3. The 

curve marked "Straightforward, low cutoff" results from perturbing Hq by V4, with operators in 
the number representation cut off at 20. The curve marked "W.r.t. local. . . " takes Hq + V2 as 
the "zero-order" Hamiltonian and perturbs it by V4 ^i^us 2 with the same cutoff, 20. This parallels 
our method for the ring problem. Finally the curve marked "large cutoff" results from perturbing 
Hq by V4 but with a cutoff of 100, yielding reliable results for levels 1 through 20. Note that the 
graph shows the (natural) logarithm of the energies (with Eq close to 1 for this A). 
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FIG. 14: Comparing spectra with and without ivi. 
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FIG. 15: Zeroth component of the local mode with and without Wi. 



Eq. ()44|) ]. and let be its annihilation operator. The frequencies, ili, are not the same as 
those defined in Eq. (fT3j) : those are for the chain, these for the ring. (When they appear 
together we distinguish using superscripts [cf. Eq. (jSD)])- We label the local mode "0" 
(although it is of highest frequency). As usual, Hq becomes 



(45) 



where we now restrict ourselves to modes symmetric under P, and call Ns the number of 
such modes. The perturbation, given earlier, is 



V 



12 2 
2^1% ) 



SO that all we need in order to proceed is the fact that 

1 



Xo 



(46) 



(47) 



Note that it is the zeroth component of each mode function, u^^\ that is important. (Had 
the asymmetric states been included they would now drop out since their zeroth components 
vanish.) In our general discussion we indicated that choosing to perturb around the local 
mode lessens the impact of the quartic term. In Eq. ^7\i the specific mechanism of that effect 
can be seen: it is built into the function Uq \ which is plotted in Fig. ^| Note that without 
uji (i.e., no local mode) the values of u'^^ fluctuate around 0.2 {u is square normalized), while 

with a moderate ui all but u^^ drop to much smaller values. 

With these deflnitions and observations the operator H can be generated. It is a sparse 
matrix and selected eigenvalues and eigenvectors can be obtained for quite large dimension. 

Our objective is to flnd the eigenstates corresponding to breathers. The most convenient 
for our purposes is the first excited state, since it presumably does not have enough excitation 
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TABLE I: Principal components of the breather state. The first 4 columns refer to the 4 symmetric 
phonons in a 6-atom ring. First row: frequencies. Subsequent rows: number operator values. 
(The highest frequency is the local mode.) Fifth column: norm squared of the mode (power of 10 
multiplier in parentheses). The local mode cutoff was 13, for the others 6. For this state first or- 
der perturbation theory is good to 0.3%. The last column reports the same calculation with cutoff 8. 
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1.728 (-5) 
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3.670 (-6) 


3.627 (-6) 



for the cutoff to be sensed (cf. Fig. lT!^ . As a state around which to perform the perturbation 
we use |0, . . . , 0, 1), which is the eigenstate of the local mode Hamiltonian, Hq of Eq. (0), 
having a single excitation of the local mode, with all other phonons in their ground state. 
The values of the parameters that we use are X = 8, uq = 1, and Ug = 1, which after 
adjusting for differences of convention, are the values used in For ui we used 2.5, 
although the structure of the eigenfunction was not sensitive to this. In Table |l] we show 
the i^o-eigenstates with the largest components in the true eigenstate. Clearly the local 
mode dominates. The next largest component is in fact the thrice excited local mode, which 
does not represent a spreading of amplitude, but rather a readjustment of the shape of 
the excitation. Other modes barely make to the 10^^ level. Of particular interest is the 
observation that the highest excitation level for other phonons is 3, indicating that a cutoff 
of 6 is safe. In fact, even to probability levels of 10~^ there is no excitation higher than 3 
for anything but the local mode. 

As a check of cutoff sensitivity we repeated this calculation with a cutoff of 8 (but with 
the local mode still at 13). The results are in the last column of Table IH There is little 
sensitivity to the change — not only in the probabilities but in the composition of the state. 

We also studied larger rings, but with smaller cutoffs. For eight atoms there were 5 
phonon modes, with a cutoff of 7 for the local-mode phonon and 5 for the others. Results 
with uJi = 2.3 are shown in Table ITTl Once again the principal change from the unperturbed 
state arises from the excited breather state. And again, down to 10~^ probability no non- 
local-mode phonon has excitation greater than one, an indication that the cutoff is again 
not felt. It is interesting that in other runs with different values of ui most of the change in 
the state was in the local-mode-phonon's contribution, which makes sense since it may be a 
better or worse approximation to the true nonlinear state, depending on Ui. 

Moving on to 12 atoms, the number of independent phonons is 7 and a cutoff of 3 was 
imposed on all but the local mode, for which the cutoff was 5. The results, in Table IIIH 
continue to be consistent with our previous conclusions. 

Of particular interest for our own application is the case of nonlinearity in the nearest 
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TABLE II: Principal components of the breather state. The first 5 columns refer to the 5 
symmetric phonons in a 8-atom ring. First row: frequencies. Subsequent rows: number operator 

values. (The highest frequency is the local mode.) The sixth column is the norm squared of the 
mode. The local mode cutoff was 7, for the others 5. For this state first order perturbation theory 
is good to 0.3%. 
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TABLE III: Principal components of the breather state. The first 7 columns refer to the 7 
symmetric phonons in a 12-atom ring. First row: frequencies. Subsequent rows: number operator 
values. (The highest frequency is the local mode.) The eighth column is the norm squared of the 
mode. The local mode cutoff was 5, for the others 3. For this state first order perturbation theory 
is good to 0.2%. 
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TABLE IV: Principal components of the breather state with nearest neighbor nonhnearity. Run 
parameters are followed by a description of the state according to the scheme of the previous 
tables. For each parameter set two states are described, the ground state and the state of a single 
breather excitation. 
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neighbor interaction. The results of this calculation support what we have already seen: the 
breather is very well approximated by phonons of the (quadratic) local mode. In Tables HVl 
and we show a small variation on the material displayed for the self interaction. Both 
the singly excited breather (more precisely, the true eigenstate closest to the singly excited 
local mode) and the ground state are shown. A variety of A and ui values are also used, to 
give a richer idea of variation of the state with changes in parameters. 

We summarize: in all cases the (quadratic) local mode dominates the true eigenstate of 
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TABLE V: Continuation of Table ITVl 
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the Hamiltonian, indicating, because that (quadratic) local mode is itself localized, that the 
stable eigenstates, or at least those generated by perturbation around the localized modes, 
are localized. Matching all these properties, including frequencies, we conclude that the 
quantized breather is stable. 

VI. TIME EVOLUTION 
A. General discussion of decay 

Although we have shown the breather eigenfunction (of the full Hamiltonian) to be dom- 
inated by the local modes, nevertheless there were small — very small — contributions from 
ordinary phonon states (in the local mode basis). Does this imply that a system initially 
in a state close to the eigenfunction (for example, in a local mode phonon eigenstate), must 
ultimately decay? An argument in favor of this perspective would be that so long as the 
initial state has matrix elements that connect it to a continuum of global phonons, it must 
have a decay rate, simply by Fermi's Golden Rule. And a decay rate implies exponential 
decay, even if the multiplier of time in that exponent is small. 

As a general guide the Golden Rule has wide apphcability, but it should not be made 
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into a false idol (a Golden Calf). Several years ago, in response to anomalies in computer 
decay studies, one of the present authors in collaboration with B. Gaveau [3| studied one 
kind of breakdown. The context is a state that has energy h and is coupled to a continuum 
with energies in the interval [Ei,E2], with Ei < h < E2 < oo. The level with energy h 
is within the interval and has nonzero coupling under the full Hamiltonian to the other 
levels. It was found that depending on the threshold behavior of the coupling, its strength, 
and the proximity of h to an edge, the system need not decay exponentially. What can 
happen is that when the interaction is turned on the initial state loses a bit of amplitude 
to other modes, but it eventually settles into an asymptotic state of norm well away from 
zero. To be precise, let the initial state (of energy h) be |0) and the full Hamiltonian be H. 
The survival probability is then p{t) = \{0\exp{—iHt/h)\0)\'^. Asymptotically, p{oo) (or a 
smoothed average) could be anywhere in (0, 1). 

In the emphasis was on the mathematical mechanism, namely the formation of what 
could be called a plasmon mode js^ . Suppose a single mode is coupled to levels, so that 
neither that mode nor the N levels are eigenfunctions of the full Hamiltonian. In the usual 
decay scenario, when that single mode is expanded in eigenstates of the full Hamiltonian 
its amplitude in each is of order I/a/ZV. Evolving in time, if you wait long enough there is 
"nothing" left in the original state; in the limit N 00 the decay is complete. But it can 
also happen, as it does for plasmons and for the systems studied in [l^, that the original 
state has significant [not 0(1/ \/N)] overlap with a true eigenstate, and as a result only 
decays to a value dependent on the value of that overlap. 

We illustrate both the gradual and catastrophic failure of the Golden Rule for the (2A^ + 
2) X (2A^ + 2) Hamiltonian 



where C is a column vector of coupling coefficients, /i G M, and fl = diag (co'_^, . . . u;^^), 
with Un G [Ei,E2] {h = 1). (Every decay system can be brought to this form.) Since 
threshold behavior will prove to be our main preoccupation, behavior at the other end of 
the energy range is irrelevant. For convenience we take the energy range to be [—E/2, E/2] 
and allow h to approach E/2 from below. The Un are taken to be uniformly spaced — what 
matters for decay (also when going beyond the Golden Rule) is the product p|7p with p 
the density of states and 7 = C/^/N, which is the appropriate scaling with A^. Therefore 
we can make p a constant [(2A^ + 1)/E], and build dimensional or other density of states 
features into C. We shall take C of the form C = (c/ \/N)(f), with c a constant and of the 
form 0„ = [1 — {n/{N + 1))^]^, with n = —N, . . . , N, and 6 another constant. If h is placed 
near level-n (with its associated u = u;„), the Golden Rule prediction for the decay rate is 
= 27rp(a;) 17(0;) p, which in this case becomes Tq^ = 27r(2/E)|c|2|0„|2 12^]. 
In Fig. El we show how the Golden rule prediction gradually declines in accuracy as 
h — s> E/2. For all points plotted, there is exponential decay to great accuracy and to values 
of p that reach 10""^ or less. (As lifetimes increase, A^ is also increased to allow the calculation 
to go to long times and avoid a Poincare recurrence. For the same c, 6 and h, changing 
A^ had little effect. More on this later.) Finally, for the parameter values stated in the 
figure, at about h = 0.92i?/2 a plasmon mode develops. At this point exponential decay 
is completely lost. In Fig. El we show decay behavior for h = 0.9E/2, side-by-side with 
the corresponding graph for h = 0.94, when a plasmon exists. With the formation of this 
excitation, a significant fraction of the initial amplitude never decays. This has nothing to 
do with Poincare recurrence. The same asymptotic probability as well as overall pattern of 




(48) 
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FIG. 16: "Measured" lifetime (from computer runs) compared to the Golden Rule prediction as 
the energy parameter h approaches the band edge (which is /i = 1 here). At /i = 0.92 a plasmon 
forms and "lifetime" loses meaning. Parameter values were E = 2, c = 0.2, 6 = 1/2, N = 200 for 
h ^ 0.85, then increasing values until N = 500 (matrix size is 2N + 2). Each run was long enough 
for p to reach ~ 10~^ (or smaller) with the logp-vs.-t plot remaining a straight line. 




FIG. 17: Comparing decay without and with a plasmon (h = 0.90 and 0.94, respectively). Note 
the different y-axes ranges. The plasmon energy is only slightly outside the band: the band edge 
is at 1 and the plasmon energy is 1.016. 



the curve obtains whether we did the calculation with dimension 102 or 1002 (A^ = 50 or 
500, and several values in between). This is the continuum behavior calculated in jl8[. 

Remark: Although the slight shift in h above creates markedly different long-term 
behavior, for short time the quantum Zeno effect for the two parameter sets is similar. In 
Fig. ^Jwe show in greater detail the "Zeno era," t < t^, with what I have called the 
Zeno time EUli, = h/y/{0\m\0) - {0\H\Oy. 

Remark: Formally the structure of the plasmon Hamiltonian and the structure of our 
Hq and Vj [Eqs. and Eqs. ((ZHHl)] are similar. The matrix for finding the classical local 

modes has the form Eq. ()48|1 if one first diagonalizes the "chain" as we do in Sec. IIVI The 
adjustable parameter sitting in the 1-1 position of the Hamiltonian of Eq. (|48|) is then uf. 

Next we check that the decay inhibition due to the plasmon is independent of the matrix 
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FIG. 18: A short-time close-up of the decays in Fig. 1171 For very short times one can see parabohc 
behavior (just barely discernable) and for longer times there is an anti-Zeno effect. All these 
peculiarities are phased out by about r^, the Zeno time, defined in the text. 
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FIG. 19: Logarithm of survival probability with a plasmon for different matrix sizes. As is evident, 
the asymptotic probability is well defined, with the only affect of matrix dimension being the time 
at which finite-size recurrence begins (the quantum analogue of the Poincare recurrence associated 
with finite-dimensional quasiperiodicity). 

size. This is shown in Fig. ^[ where the only change with increasing is the deferral of 
the Poincare recurrence. 



■ 102 

• 202 

• 1002 




B. Quantum time evolution of the breather 

We turn to the time dependence under quantum evolution of the breather state. Specifi- 
cally, we study \{ipo\ Gxp{—iHt / h)\ilJo)\'^ for ipo a local mode phonon (in particular, |0, . . . , 1)) 
and H the full Hamiltonian. This is shown in Fig. 1201 |i3^] . As we saw in the examples of 
Sec. IVl Al where a plasmon has formed, initial decay is followed by stabilization bounded 
away from zero. In the present situation the value at which it stabilizes is much larger than 
the ~0.6 (cf. Fig. ITIJ|) of those examples, which reflects the much larger amplitude of the 
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FIG. 20: Survival as a function of time; the initial state is the singly excited local mode. 

local mode eigenstate in the true eigenstate of the Hamiltonian (cf. Table H]). Although the 
coupling pattern for the two cases is not the same (we will turn to this in a bit), the basic 
idea is the same: the initial state has order unity overlap with a true eigenfunction, and 
although it has coupling to the continuum, that coupling does not lead to exponential decay. 



Since for the plasmon mode the emphasis is on the threshold structure of the coupling, 
it is of interest to explore those properties for the breather. This is also useful for assessing 
properties of the eigenstate when, because of the size of the system, numerical diagonalization 
is out of reach. 

States of the local mode are both our initial states for decay and our unperturbed states 
for the numerical diagonalization. One of these is |local) = |0, 0, 0, . . . , 1) where, as in Sec. 
|VJ the "1" refers to the eigenvalue of the local mode number operator and the zeros to other 
phonon levels. The couphng, the analogue of "C" of Eq. (jiH|) . is 



In Fig. 1^ is a logarithmic plot of (? as a function of energy. With states ordered by energy 
[as in Eq. (j48|) ] this clearly does not provide a smooth function. However, the principal 
issue is not a precise resemblance to the plasmon paradigm. The plasmon was discussed 
not because of its specific coupling pattern, but rather as an example of how, although the 
initial state does couple to a continuum, the initial state still has finite (i.e., not going to 
zero with system size) overlap with a true eigenstate of the system. The pattern seen in 
Fig. 1^ shown for 14 oscillators, is the same as one sees for 4, 6, 8 or any number that we 
have been able to study. As such, the substantial overlap of the true eigenstate (induced by 
the nonlinear A) with the local mode (induced by ui) should continue as system size grows, 
since it these matrix elements that determine the overlap. 

Remark: Note that stability may not persist in all dimensions; certainly threshold 
features of the density of states and spectrum are affected by dimension, and the usual 
intuitions regarding Fermi's Golden rule may again hold sway. In 2] we made the point 
that the symmetry breaking of the Jahn- Teller effect makes this a one-dimensional problem. 



1. Off-diagonal matrix elements of the singly- excited local mode 



g{\ni, ^2, . . .)) = (local|a;o|ni, ns, . . .) . 



(49) 
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FIG. 21: Base-10 logarithm of | (0, 0, . . . , 0, l|xQ|ni, n2, . . .)| for a ring of 14 oscillators. States, 
|ni, 71-2, . . .), are ordered by increasing energy. 

significantly enhancing the possibility of classical breathers. The same may be true quantum 
mechanically. It may even be that this plays a role in the temperature-dependent decay 
of the breather (through an effective increase in dimension), as evidenced by the high- 
temperature disappearance of anomalous decay in doped alkali halides i^]. For numerical 



diagonalization, increase in dimension is difficult computationally. However, for the path 
integral there should be little problem, and indeed 3-dimensional phonon coupling is used 
in many significant problems (see for example [30,]). 

VII. CONCLUSIONS 

Both numerical diagonalization and path integration automatically deal with infinite- 
lifetime eigenstates of the Hamiltonian. The principal issue is therefore to show localization 
for the system we study. That system is not the full periodic system, but rather a version 
that neglects quantum-tunneling delocalization [an order exp(— [positive constant]//i) effect], 
but otherwise differs little from the full translationally invariant system. 

Localization is demonstrated in two ways. For the path integral we show that increasing 
value of the quartic coupling constant ("A") decouples the rapidly oscillating atom from 
distant atoms, in the same way the this phenomenon occurs for a quartic localized mode. 
For the method of numerical diagonalization we show that the breather state is, up to tiny 
corrections, entirely constructed of states of an appropriate quadratic local mode. We remark 
that there is no principle saying that such tiny corrections represent either delocalization 
or decay. For example, the harmonic oscillator state, iIj{x) ~ exp(— mu;a;^/2) is certainly a 
localized object, even though it is nonzero for all x. As to decay, we devote an entire section 
to disabusing anyone of the notion that decay is a mathematical imperative, although what 
lies behind the non-decay is merely the fact that under the full Hamiltonian the overlap of 
a state ipo with itself does not go to zero if ipQ has finite (bounded away from zero) overlap 
with an eigenstate of the full Hamiltonian. 

It follows that the quantized breather is stable, except against spontaneous tunneling, 
a process that is in general far slower than the decay rates that have been ascribed to the 
breather. 
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